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SECTION  I 


DC  RESTORATION 


INTRODUCTION 

This  report  is  the  first  interim  progress  report  for  Contract  DAAG53-76-C- 
0195,  Automated  Image  Enhancement  Techniques  for  2nd  Generation  FLIR. 

It  describes  results  of  Phase  I studies  of  image  processing  applications  to 
FLIR  imaging  systems  from  July  1,  1976  to  December  1,  1976. 

The  report  is  subdivided  into  five  major  headings  corresponding  to  Phase  I 
program  tasks: 

• DC  Restoration 

• Contrast  Enhancement 

• MRT  Enhancement 

• Resolution  Restoration 

• Statistical  Analysis 

DC  restoration  treats  the  problems  of  correcting  or  compensating  for  defic- 
iencies in  detectors  and  signal  processing  electronics  which  can  result  in 
harmful  distortion  of  the  IR  image.  Contrast  enhancement  deals  with  prob- 
lems of  transforming  thermal  images  of  wide  dynamic  range  into  the  limited 
dynamic  range  of  the  display  medium  without  losing  essential  information 
content.  MRT  enhancement  involves  spatial  or  temporal  smoothing  for  noise 


reduction.  Resolution  restoration  treats  detail  recovery  from  imagery  which 
has  been  blurred  by  the  optics  and  electronics  of  the  sensor.  Lastly,  statis- 
tical analysis  of  FLIR  imagery  is  performed  to  provide  quantifiable  charac- 
terizations through  shape,  brightness  and  textural  feature  statistics. 

I 

• I 

DC  RESTORATION 


Direct  versus  AC  Coupling 

DC  restoration  traditionally  refers  to  techniques  for  correcting  imagery 
artifacts  introduced  by  AC  coupling  detectors  to  their  amplifying  electronics. 
The  reasons  for  AC  coupling  (usually  with  simple  high  pass  RC  networks) 
are: 

• Background  suppression  for  contrast  rendition  and  efficient 
utilization  of  electronic  dynamic  range 

• Removal  of  detector  DC  bias 

• l/f  noise  suppression 

Artifacts  created  by  AC  coupling  are  transient  effects  (drooping  across  the 
scan  line  and  undershooting  on  hot-cold  transitions)  and  information  loss 
about  inter  channel  DC  levels  occurs.  All  channels  have  zero  average 
signal. 

The  most  common  DC  restore  technique,  which  corrects  the  second  defect, 
clamps  the  signals  on  all  channels  by  imaging  a common  reference  source 
prior  to  each  scan.  Each  channel  then  measures  scene  radiance  changes 
relative  to  the  same  reference  level,  and  relative  interchannel  DC  levels 
are  preserved. 


8 


With  this  approach  it  is  difficult  to  adapt  the  reference-source  temperature 
to  the  ambient-scene  temperature.  As  a result,  the  dynamic  range  of  the 
electronics  may  be  underutilized,  and  needless  saturation  at  the  hot  or 
cold  end  of  the  range  will  occur  if  the  respective  references  are  too  hot  or 
too  cold. 

With  the  advent  of  monolithic  and  hybrid  focal  plane  detector/CCD  proces- 
sors, it  is  uncertain  whether  or  not  the  previous  reasons  for  AC  coupling 
detector  outputs  are  still  valid. 

Directly  coupled  PV  IR/CCDs  exhibit  negligible  1/f  noise  (discussed  later 
in  this  section).  This  was  one  of  the  important  reasons  for  AC  coupling 
when  using  PC  detectors.  Since  a compromise  between  the  deleterious 
effects  of  1/f  noise  and  AC  coupling  transients  are  no  longer  needed,  RC 
time  constants  that  are  long  relative  to  a line  scan  interval  (RC  > 1 sec  for 
100  channels  and  < 1%  droop  across  a scan  line)  are  now  desired.  On  chip 
capacitors  up  to  about  100  pf  are  practical,  but  fabrication  is  difficult 
when  in  series  rather  than  shunt.  Also,  RC  > 1 sec  required  R >10^^  ohms, 
which  may  be  difficult  to  achieve  as  the  input  impedance  to  the  CCD  input 
amplifier. 

Adaptive  background  subtraction  ideally  ensures  that  the  coldest  part  of  the 
scene  is  always  at,  or  slightly  above,  the  zero  signal  level.  Prevailing 
direct  coupled  PV  IR/CCD  techniques  handle  this  at  the  CCD  input  stage  by 
variable  gate  thresholding  which  transfers  only  the  charge  above  a background 
level  determined  by  the  previous  frame.  A comparable  adaptive  approach 
for  a clamped,  AC-coupled  system  would  require  injecting  a scene  dependent 
charge  directly  on  the  coupling  capacitor  during  the  clamping  interval. 
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rather  than  imaging  a constant  reference  temperature.  By  so  doing,  however, 
the  advantage  of  simple  cancellation  of  detector  DC  bias  through  AC  coupling 
would  be  lost.  Rather,  as  in  the  DC-coupling  case,  the  detector  bias  for 
each  channel  would  have  to  be  added  to  the  background  voltage  to  provide 
both  corrections  simultaneously. 


Because  the  focal  plane  CCDs  are  cooled,  a dynamic  range  of  80  dB  (10,  000:1 
peak  signal/rms  noise)  is  expected.  In  order  to  estimate  the  scene  temper- 
ature range  which  can  be  accommodated  by  the  electronics  on  the  focal  plane, 
we  compute  the  spectral  radiance  in  the  8-12  /um  band  (from  Planck's 
radiation  formula).  This  indicates  that  a spectral  radiance  resolution  of 
5R  = 0.022  watts/m  sterrad  is  equivalent  to  a 0.  05°K  temperature  resolu- 
tion at  a nominal  ambient  temperature.  At  a minimum  temperature  of 

2 

T . = 200°K;  for  example,  R . = 4W/m  sterrad.  Then  for 

min  min 

R = R + 10^6R  = 224,  T = 475°K. 
max  min  max 


Assuming  that  all  gain  factors  (detector  size,  IFOV,  detector  responsivity, 
input  amplifier  gain,  etc. ) are  such  that  R translates  into  the  maximum 
CCD  well  charge,  the  focal  plane  electronics  is  capable  of  handling  an 
equivalent  scene  temperature  range  of  200°-475°K  (-100°F  to  395°F)  at 
0.05°K  resolution.  Raising  the  minimum  temperature  to  293°K  only  raises 
the  maximum  temperature  to  490°K  (68°F-420°F). 


This  calculation  illustrates  two  important  points  with  regard  to  background 
subtraction.  First,  the  dynamic  range  of  cooled  CCDs  is  probably  sufficient 
to  permit  a constant  background  level  subtraction  and  still  preserve  all 
significant  information  in  a typical  scene.  Second,  the  nonlinear  relation- 
ship between  temperature  and  observable  radiant  power  quickly  produces 
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I 

diminishing  returns  in  equivalent  temperature  dynamic  range  (at  fixed  reso- 
^ lution)  as  the  background  subtraction  level  is  raised.  The  off-focal  plane 

electronics  outside  the  dewar  will  have  dynamic  range  of  about  60  dR,  and  it 
will  lead  to  a display  with  a somewhat  less  dynamic  range.  Thus,  it  appears 
much  more  important  to  provide  on-focal  plane  contrast  enhancement  to 
ensure  that  the  information  a /ailable  inside  the  dewar  is  preserved  when  it 
I is  taken  outside. 

t 

I 

In  summary,  a number  of  PV  IR/CCD  fabrication  and  performance  issues 
need  to  be  resolved  before  answering  the  question  of  whether  or  not  AC 
coupling  on  the  focal  plane  should  be  implemented.  Examination  of  these 
issues  by  detector  and  focal  plane  electronics  designers  at  IIRC  has  begun. 

In  any  event,  it  is  probably  desirable  to  incorporate  AC  coupled  electronics 
off  the  focal  plane  to  avoid  the  problems  of  drift  in  DC  coupled  amplifiers. 
Here,  long  time  constant  coupling  networks  are  easily  realized,  and  channel 
DC  levels  established  on  the  focal  plane  can  be  preserved  simply  by  clamping 
to  zero  at  the  start  of  every  scan  line. 


Uniformity  of  channel-to-channel  responsivity  is  an  important  requirement 
if  high  temperature  resolution  is  to  have  any  meaning.  PV  detectors  can  be 
fabricated  at  about  5%  uniformity,  but  temperature  resolution  of  order  0.05°K 
require  uniformity  of  0.  1%  or  better  among  channels. 

In  parallel  scanned  systems  with  serial  TDI,  there  is  a question  about  the 
difficulty  of  providing  a spatially  uniform  calibration  source  for  the  detector 
array.  In  this  section  we  examine  a signal  processing  alternative  that  does 
not  rely  upon  reference  source  uniformity. 
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A TDI  channel  can  be  treated  as  a single  detector  provided  that  each  detec- 
tor in  the  channel  TDI  array  views  the  same  spot  on  the  reference  (as  it 
does  during  scanning  the  scene).  Thus,  if  is  the  responsivity  of  the  j 
detector  in  the  channel.  is  its  DC  bias,  and  P^.(t)  is  the  radiant 
power  at  time  t.  the  output  voltage  of  the  detector  is: 


kj 


The  channel  output  for  (J+1)  TDI  detectors  is 


If  each  detector  views  the  same  point  in  sequence  at  intervals  6, 

- (J-i)6]  = Pj^j(t)  = Pk(t).  0 < j < J 

Hence,  we  can  write 


' '’kPk'*>  ®k 


where 


R. 


J 

I 


j=0 


J 

z 


j=0 


Channel  responsivity  equalization  can  therefore  be  treated  as  single  detec- 
tor equalization. 
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The  calibration  technique  assumes  that  a pair  of  line  sources  at  different 
temperatures  (T,.  T„),  conveniently  located  on  the  field  stop  for  example, 

± A 

are  sensed  at  the  end  of  each  line  scan.  The  difference  signal, 

W - '^2'  - vv^i” " ^ w * w 

is  taken  to  remove  the  bias  (Nj^  is  the  noise).  At  successive  field  intervals, 
T,  the  vertical  scan  mirror  is  nodded  one  IFOV  so  that  adjacent  detectors 
have  imaged  the  same  reference  spot  as  illustrated  in  Figure  1. 


Reference  Channel 


t.  + t: 
1 


•''l 

-K^.l 

-2 

-1 

0 

1 

-1 

0 

t 

2 

K„ 


Figure  1.  Overscan  Geometry 


Calibration  of  the  array  is  performed  in  pairs --channel  1 relative  to  channel 
0,  2 relative  to  1,  etc.  All  channels  are  referenced  to  the  center  channel 
since  errors  propagate  in  direct  proportion  to  distance  from  the  reference 
channel.  A simple  iterative  algorithm  of  the  form  (for  k>0) 


^^-2 
B=  c/  L zf 

k=0  ^ 
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permits  continuous  on-line  adaptation  to  changing  environment  of  both  the 
detectors  and  calibration  source.  Although  a direct  adaptation  of  the  respon- 
sivity  is  indicated,  in  actuality  the  gain  of  the  channel  output  amplifiers 
would  be  controlled  to  ensure  equal  gains  for  the  full  channel  electronics. 

Figure  2 shows  results  of  a simulation  of  adaptive  gain  equalization  for  15 
elements.  The  logj^  of  the  error  in  equalized  gain  for  the  first  detector 
(1-dots)  and  the  center  detector  (8-+  signs)  and  the  last  detector  (15-solid 
curve)  relative  to  the  reference  detector  (0)  are  all  plotted  for  each  of  100 
iterations  (1.  7 sec  at  a field  update  rate).  The  S/N  is  30:1  (real  ratio). 

The  initial  detector  res  pons  ivities  for  all  detectors  except  1,  8,  and  15  were 
determined  by  a random  number  draw  from  a Gaussian  distribution  (+  10% 
one  sigma);  and  the  calibration  source  was  assumed  to  have  5%  one  sigma 
variability  from  spot-to-spot  and  look-to-look.  Initial  gain  errors  of  +30% 
for  detector  1,  -30%  for  detector  8,  and  +20%  for  detector  15  were  chosen 
to  test  the  convergence  of  the  algorithm.  The  iteration  parameter  constant, 
c,  was  decreased  linearly  from  0.  3 down  to  0.  03  after  50  iterations  to  speed 
convergence,  thereafter  it  was  held  constant  at  0.03.  After  100  iterations 

logjQ|(R^/R^)  - 1|  <-3,  l<k<  15 

indicating  that  all  channel  gains  were  equalized  to  within  0. 1%. 

Experience  has  shown  that  this  scheme  fails  for  very  long  arrays.  How- 
ever, by  dividing  the  array  into  subarrays,  the  maximum  calibration  dis- 
tance to  the  last  detector  in  any  subarray  can  be  minimized.  A scheme  for 
equalizing  the  maximum  calibration  distance  is  illustrated  in  Figure  3 for 

the  case  N=10,  K=4.  The  lead  detector  in  each  subarray  (k  ) is  calibrated 

o 
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Figure  2.  Adaptive  Responsivity  Equalization 
Simulation  Results 

with  respect  to  the  reference  detector  (0  ),  then  the  other  subarray  detec- 
tors  are  calibrated  with  respect  to  their  lead  detector.  This  approach 
requires  K overscans  to  obtain  all  the  information  required  for  each  iter- 
ation. In  Figure  3,  the  first  overscan  provides  the  single  detector  shift 
required  for  calibration  within  the  subarrays,  the  remaining  provide  the 
overscan  of  the  subarray  references;  i.  e. , 3^  with  2^  at  t.  + 2t,  2^  with  1^ 

at  t + 3t,  and  finally  1 with  0 at  t.  + 4t.  Although  this  reduces  the  iter- 
1 o o 1 

ation  rate,  different  overscans  could  be  performed  on  successive  line  scans, 
rather  than  successive  fields,  if  less  than  full  parallel  scan  is  implemented. 
Note  that  for  100  channels,  N=50  and  K=10.  Thus,  no  channel  would  be 
more  than  9 away  from  the  reference,  and  rapid  convergence  would  be 
expected  based  upon  the  results  for  K*=15  in  Figure  2. 
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Figures,  Equal  Calibration  Distance  Subarrays:  N<K(K+l)/2 

As  a final  comment,  PV  detectors  exhibit  good  linearity  between  radiant 
power  and  detector  current  over  several  decades  of  incident  radiance.  This 
holds  true  for  detectors  coupled  to  CCDs  as  well,  provided  that  background 
current  bias  offset  is  incorporated  to  maintain  a detector  resistance, 

CCD  input  amplifier  transconductance  product  greater  than  one 
The  dynamic  range  example  discussed  previously^  indicated  a radiance  range 
of  only  about  55:1  for  full  80  dB  CCD  dynamic  range  utilization  at  0.05°K 
resolution.  Therefore,  it  appears  that  the  calibration  reference  tempera- 
ture need  not  be  adapted  to  the  ambient-scene  temperature  to  provide  a good 
calibration  over  the  full  dynamic  range  of  interest. 
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PV  IR/CCD  l/f  Noise 


It  was  stated  previously  that  l/f  noise  in  direct-coupled  PV  IR/CCDs  was 
negligible.  This  section  presents  some  recent  experimental  evidence  typical 
of  hybrid  (FbSnTe  and  HgCdTe)  detector /CCDs  which  supports  that  state- 
ment. Figure  4 is  a plot  of  experimental  noise  spectral  data  at  60T|A  CCD 
input  amplifier  drain  current.  Here  again,  a high  product  produces 

low  drain  currents  and  correspondingly  low  l/f  noise  rules.  Actually,  30T]A 
for  8-12um  and  3T\A  for  3-5um  are  more  typical. 


Figure  4.  PV  Detector/ CCD  Noise 


For  a noise  spectrum  of  the  form 

P (f)  = C(f.  /f+1) 
n k 

in  a passband  of  f,  < f < f„,  the  noise  variance  associated  with  frequencies 

below  f,  is 
k 

and  with  frequencies  above  f^^ 

2 

We  know  that  we  cannot  take  fj=0  in  this  model  since  = “ would  result. 

Since  l/f  noise  has  been  observed  down  to  10  ^ Hz  (2.  8 hour  time  constants), 
assume  f^  = 10"*^  Hz  is  a lower  bound  beyond  which  the  low  frequency  noise 
essentially  has  constant  spectral  density,  and  hence  negligible  contribution 
to  the  total  energy. 

In  an  875  line  TV  compatible  system  with  2:1  interface,  the  field  scan  time 

at  83%  scan  efficiency  is  about  13.  8 msec.  A 100  channel  system  requires 

an  information  bandwidth  of  about  fg  = 184  kHz  for  4:3  aspect.  At  f^  - 10  Hz, 

f,  = 250  Hz  and  f„  = 184  kHz, 
k 2 

Oj  = (3.9  kHz)C 
a|  = (185.4  kHz)C 
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/ 2 2 

Since  the  total  noise  standard  deviationTyaj  + <^2  ' about  1%  higher 

than  a„,  it  is  extremely  unlikely  that  any  low  frequency  flicker  noise  would 

be  noticeable.  In  the  worst  case  of  a fully  parallel  scan  system,  f.,‘^42  kHz 

#2  2 ^ 
andyaj  + is  about  4.  5%  higher  than  indicating  l/f  noise  would  still 

be  neglibible  at  the  slowest  scan  rate. 


SECTION  II 


CONTRAST  ENHANCEMENT 


INTRODUCTION 

This  task  addresses  two  related  problems  in  thermal  imaging: 

1.  The  dynamic  range  of  the  temperatures  within  a scene 
(typically  1000:1)  can  be  far  higher  than  the  display  dynamic 
range  (less  than  100:1). 

2.  Moreover,  the  scene  dynamic  range  itself  can  change  from 
scene  to  scene,  requiring  the  global  gain  and  brightness 
controls  to  be  scene  adaptive. 

The  limited  display  dynamic  range  results  in  the  loss  of  local  contrast 
masking  target  details.  The  objective  is  to  reduce  the  dynamic  range  over 
the  scene,  while  maintaining  (or  even  improving)  the  local  contrast  in  the 
image.  Further,  the  scene -to -scene  brightness  range  variation  should  be 
compensated  for  by  an  automatic  global  gain  brightness  control.  In  this 
section,  various  local  contrast  enhancement  and  automatic  gain  brightness 
control  techniques  being  applied  to  the  FLIR  imagery  are  reported. 

We  assume  the  display  characteristics  to  be  linear  between  the  pedestal 
brightness  B and  the  dynamic  range  given  by  B + R.  Given  a scene  of 
arbitrary  dynamic  range  Rg  and  the  pedestal  brightness  Bg,  we  have  to 
transform  the  input  image  to  fit  the  range  B and  B + R.  Since  Rg  is  usually 
much  larger  than  R,  this  involves  compression  of  the  dynamic  range  of  the 


scene.  Local  contrast  would  be  lost  if  linear  compression  were  used. 

Hence,  we  must  make  better  use  of  the  available  grey  levels  in  the  display. 
Second,  the  relative  brightness  of  distant  areas  within  the  scene  contributes 
little  to  discrimination.  Local  area  brightness  controls  can  help  by  sup- 
pressing slow  variations  of  the  brightness  across  the  image,  while  retaining 
the  local  contrast.  This  concept  immediately  leads  to  high-emphasis  linear 
and  homomorphic  filteiing,  successfully  employed  for  image  enhancement. 
Indeed,  we  will  establish  the  correspondence  between  the  two  approaches  in 
a later  section. 

The  techniques  being  investigated  are; 

• Local  area  histogram  modification 

• Local  area  gain  brightness  control  adaptive  to  global  range 
variations 

• Linear  and  homomorphic  filtering 

These  techniques  are  not  mutually  exclusive.  The  goal  of  this  study  is  to 
define  a scheme  that  combines  the  best  features  of  each  approach  applicable 
in  the  FLIR  context. 

HISTOGRAM  MODIFICATION 

Histogram  modification  refers  to  a non-linear  mapping  of  the  original  inten- 
sities to  another  domain,  so  that  the  new  intensities  have  a specified  distri- 
bution. When  the  new  distribution  is  uniform,  the  mapping  process  is 
referred  to  as  histogram  equalization.  The  result  of  histogram  equaliza- 
tion is  that  more  frequent  grey  levels  are  spread  farther  so  that  they  occupy 


a greater  portion  of  the  available  range,  while  the  less  populated  grey  levels 
are  drawn  closer. 


The  histogram  modification  routines  are  based  on  the  following  analysis: 


Let  the  original  intensity  levels  be  denoted  by  x = K^,  ....  K^.  The  corres- 
ponding histogram  normalized  into  a density  function  (by  dividing  by  the 
number  of  elements  in  the  area),  is  given  by  pj^(x),  x = K^,  . . . , Kg.  We 
desire  a mapping  y = g(x)  such  that  y has  a given  density  function  Py^y). 
y = N, N„. 


Let  the  corresponding  cumulative  distributions  be  P„(x)  = Ep  (x) 

2\  X 


X = K. 


and  P^(y)  = EPy<y) 

y = N, 


Since  both  Pj^(x)  and  Pyly)  are  monotonically  increasing,  it  can  be  shown 


y = g(x)  = P'^  CPx<x>] 


Figure  5 illustrates  this  mapping  schematically.  A grey  level  such  as  x^ 
is  mapped  into  new  intensity  y^  by  the  above  transformation.  Pj^(x)  then 
is  the  cumulative  distribution  computed  from  the  desired  histogram  shape. 
In  Figure  5 we  show  P (y)  to  be  a linear  ramp  corresponding  to  a uniform 
density  function.  Actually,  if  the  desired  density  is  uniform,  it  is  only 
necessary  to  look  up  the  value  Pj^(x)  for  a given  x.  The  new  intensity 
y = g(x)  is  then  simply  linearly  proportional  to  P^(x),  For  an  arbitrary 


distribution  Pyiy)  however,  the  second  look-up  implied  in  Figure  5 has  to 
be  performed  to  find  the  new  transformed  intensity  y^.  Figure  5 also  shows 
the  way  for  an  efficient  look-up  table  implementation  of  histogram  modifi- 
cation using  CCDs. 


Figure  5.  Histogram  Modification  Mapping 


Full  Frame  Histogram  Modification 

Full  frame  histogram  modifications  with  uniform  (equalization)  and  hyperbolic 
distributions  have  been  tried.  Figure  6 is  an  8-bit  (256  levels)  quantized 
NVL-supplied  image  reproduced  on  an  Optronics  film  writer  constrained  to 
a 50:1  dynamic  range.  Figure  7 is  the  corresponding  histogram.  Note  that 
there  are  only  163  levels  in  the  range.  This  image  was  then  histogram 
equalized  so  that  the  resultant  histogram  would  be  uniform  in  the  range 
0-255.  Figure  8 is  the  corresponding  equalized  image. 
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An  alternative  to  equalization  is  histogram  hyperbolization  (reference  1). 
This  is  based  on  the  premise  that  the  perceived  brightness  b is  proportional 
to  the  logarithm  of  the  display  intensity  I.  Therefore,  the  observed  bright- 
ness levels  (log  I)  should  be  uniformly  distributed.  To  achieve  this,  the 
distribution  of  the  displayed  intensities  should  be  hyperbolic,  instead  of 
uniform.  Figure  9 is  a reproduction  of  the  image  modified  in  this  manner. 
Histogram  hyperbolization  appears  to  be  an  improvement  over  equalization. 

Local  Area  Histogram  Modification 


In  local  area  histogram  modification  the  image  points  in  a small  area 
(window)  are  modified  using  a transformation  computed  from  the  histogram 
of  the  larger  subregion  surrounding  the  local  area.  Figure  10  illustrates 
the  inner  and  outer  regions.  Local  area  histogram  equalization  might  alle- 


Figure  8.  h’ull  Frame  Histogram  Equalized  Image 
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Figure  9.  F’ull  Frame  Histogram  Hyperbolized  Image 


Compute 

Histogram 


Moving  Window 


Illustrates  Inner  and  Outer  Windows  for 
Local  Area  Histogram  Modification 


viate  some  problems  of  full  frame  equalization,  because  the  grey  scale  trans- 
formation is  now  based  on  the  area  immediate  to  the  area  being  equalized. 
However,  the  following  questions  arise  in  local  area  histogram  equalization: 

The  limits  of  the  equalized  histogram  (N^^,  N2)  should  be  determined  for  each 
local  area.  If  all  local  area  histograms  are  equalized  out  to  the  same  inten- 
sity limits  as  in  Figure  11,  the  relative  contrast  between  two  areas  with 
totally  different  brightness  is  completely  lost  and  undesirable  texture  is 
brought  out  in  bold  relief  as  shown  in  Figure  12,  Indeed,  this  approach 
corresponds  to  completely  removing  the  low  spatial  frequency  components 
in  the  image--clearly  an  unwarranted  phenomenon. 

An  alternate  approach  is  to  make  the  new  limits  (N^,  Ng)  directly  propor- 
tional (Figure  lib)  to  the  old  limits  (K^^,  Kg)  in  each  local  area  as  follows: 

= aK^  + P and  Ng  = aKg  +•  p 

where  a and  p are  chosen  so  that  the  maximum  and  minimum  transformed 
intensities  (over  all  local  areas)  correspond  to  the  extreme  limits  of  the 
display  - (B,  B + R).  Figure  13  was  equalized  in  this  manner,  and  we  see 
that  the  relative  contrasts  between  areas  of  different  brightness  are  pre- 
served. However,  note  the  boxy  anomaly  around  the  hot  target,  where  the 
severe  banding  effect  has  been  accentuated.  This  is  because  the  local  area 
histogram  is  strongly  bimodal  on  a target /background  boundary  as  in 
Figure  14,  The  higher  intensity  peak  corresponds  to  the  hotter  target  and 
the  lower,  broader  peak  to  the  background.  The  transformed  intensities 
are  shown  below  the  histogram.  We  note  that  the  equalization  compresses 
the  valley  region  between  the  two  peaks  in  the  histogram.  This  raises  the 
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a.  Same  Interval  (N N2)  for  all  Local  Areas 


b.  Proportional  Intervals 

= oK^  + e.  Ng  = aK2  + P 


nan  n = A(K  - K ) + Bias 

r g • 

LAGBC  Incorporated  into  the  Equalization  Intervals 

Figure  11.  Different  Equalization  Interval  Schemes 
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Figure  12.  Local  Area  Histogram  Equalized,  as  in  Figure  Ha 


Figure  13,  Local  Area  Histogram  Equalized,  as  in  Figure  lib 


Transformed  Grey  Level  Scale 


Figure  14.  Shows  Bimodal  Distribution  on  a Local 
Area  with  Target  and  Background 

intensity  level  of  middle  intensity  background  region  (where  banding  has 

occurred)  closer  to  the  hot  target,  accenting  the  banding  effect.  To  alleviate 

this,  other  shapes,  including  the  bimodal  catenary,  are  being  investigated 

N - N 

for  histogram  modification.  However,  the  gain  for  each  local  area  1 2 

is  not  spatial  frequency  dependent  as  it  was  with  Figure  12.  We  1 2 

would  still  like  to  attenuate  low  spatial  frequencies  in  the  image  (without 
removing  them  altogether).  This  would  make  more  room  in  the  display 
dynamic  range  for  local  contrast  to  be  enhanced. 

A third  approach  is  to  make  the  ranges  (N^,  Ng)  adaptive  over  each  local 
area  in  a manner  precisely  analogous  to  the  local  area  gain  brightness 
approach  to  be  outlined  in  the  next  subsection.  We  will  defer  the  discussion 
of  this  equalization  until  then. 
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lcx:al  area  gain  brightness  control 


Given  a scene  that  has  a dynamic  range  between  Bg  and  Bg  + Rg,  local  area 
gain  brightness  control  should: 

• Vary  local  average  brightness  to  contract  the  overall  dynamic 
range  of  scene 

• Enhance  local  variations  above  the  human  contrast  sensitivity 
threshold  to  bring  out  the  details 

• Automatically  fit  the  new  extremes  in  the  image  to  the  display 
limits  B and  B + R 

The  above  assumes  that  the  dynamic  range  in  a scene  is  larger  than  that  of 
the  display.  It  is  of  course  possible  to  have  a low  dynamic  range  image 
where  it  is  not  necessary  to  compress  the  overall  dynamic  range  to  fit  the 
display.  Therefore,  this  feature  should  be  scene  adaptive.  Following  is  a 
formulation  that  does  precisely  this: 

The  image  intensity  I.,  of  each  point  is  transformed  based  on  local  area 
statistics — the  local  mean  m..  and  standard  deviation  a.,  computed  on  a 
local  area  surrounding  the  point.  The  new  intensity, 

A 

I..  = G..(I..  - m )+A(m..  -m)+H+B 
1]  ij  U H H 

V J V„  ..  --  . ^ 

(1)  (2) 

and  the  local  gain, 

G..  = [A(m..  - m)  + H + B] 

IJ  QfO..  ij 
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M is  the  global  mean  and  H is  a bias  value  computed  to  restore  the  average 
value  to  the  image.  A is  constant  for  a given  scene.  A and  II  are  determined 
by  the  actual  dynamic  range  present  in  the  image. 

This  equation,  based  upon  the  model  shown  in  Figure  15a,  shows  that  the 
transformed  image  is  composed  of  a slowly  varying  brightness  component 
(2)  and  superimposed  higher  frequency  local  variations  (1).  The  local  gain 
operates  on  the  first  term  and  is  itself  proportional  to  the  average  local  dis- 
play brightness  component  (2).  As  seen  in  Figure  15b,  the  gain  A reduces 
the  slowly  varying  component  bringing  it  closer  to  the  global  mean,  so  that 
there  is  more  room  for  enhancing  local  details. 

The  rationale  for  making  the  local  gains  G_  proportional  to  the  local  display 

brightness  is  that  over  a large  range  of  intensities,  the  minimum  resolvable 

intensity  difference  Al  is  directly  proportional  to  the  average  local  intensity 

I.  The  contrast  sensitivity  C is  the  proportionality  constant.  Also  the 

s 

local  gain  is  made  inversely  proportional  to  the  local  standard-deviation 
effort  to  enhance  areas  of  little  contrast. 

Note  that  the  above  formulation  is  similar  to  that  of  reference  2,  but  with 
the  following  constraints. 

• The  gain  A and  the  bias  H are  computed  on  a per  scene  basis 
to  fit  the  modified  scene  range  to  the  display.  A can  there- 
fore be  greater  than  unity  (for  scenes  with  low  dynamic  range) 

• The  local  gains  are  now  directly  related  to  the  contrast 
sensitivity  of  the  human  eye 

• Standard  deviation  determines  the  extent  to  which  the  local 
standard  deviation  governs  the  local  gain 
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The  size  of  the  local  area  over  which  the  local  means  are  determined  is 
important.  This  determines  the  spatial  frequencies  that  are  attenuated 
during  the  LAGBC  process.  An  exact  correspondence  between  this  and 
spatial  frequence  filtering  will  be  established  in  a later  section. 

Implementation 

The  primary  computation  involved  is  in  evaluating  the  local  means  and  stan- 
dard deviation  over  a sliding  window  centered  at  each  point.  This  is  ineffi- 
cient from  the  implementation  viewpoint.  Two  alternate  approaches  are 
being  implemented; 

1.  Similar  to  the  overlapping  window  approach  used  in  Local 
Area  Histogram  Equalization  (Figure  10).  The  mean  and 
standard  deviation  coniputed  for  the  outer  window  are  assumed 
constant  over  the  inner  window.  We  call  this  the  piecewise 
constant  approach. 

2.  The  second  approach  is  to  have  nonoverlappin;^  local  regions. 

The  local  mean  and  standard  deviation  at  each  point  are 
obtained  by  a simple  two-dimensional  interpolation.  This  is 
a piecewise  linear  approximation. 

Local  Area  Gain  and  Brightness  Applied  to  Histogram  Equalization 


Figure  11c  illustrates  this  approach  where  the  limits  of  equalization  for  each 
area  histogram  are  determined  by  considerations  similar  to  that  of  LAGBC. 


By  comparison  with  the  LAGBC  equation,  we  see  that, 

• Mean  local  display  intensity  n 

""2 

• Local  gain  is  approximated  by  ^ factor 

' 2 


The  local  gain  should  be  proportional  to  the  local  mean  brightness;  i,e. , 


"l  ■ "2 


K. 


or . 


Local 


If  we  approximate  ” ^2^’  spread  of  the  LAH,  we  get 


which  implies  that  the  range  of  the  mapped  intensities  should  be  directly 
proportional  to  the  average  mapped  intensity  for  that  local  area.  Also,  by 
comparison  with  the  LAGBC  equation,  the  average  local  display  brightness  is 

n = A(K  - Kg)  + H + B 


where  Kg  is  the  global  intensity  mean  and  B is  the  pedestal  brightness  as 
before.  A is  the  brightness  gain  and  H is  a bias,  computed  as  before,  so 
that  the  global  extrema  of  the  mapped  intensities  correspond  to  the  limits 
of  the  display  (B,  B+R).  This  approach  is  being  implemented  and  prelim- 
inary experimentation  with  different  values  of  the  proportionality  constant 
is  being  made. 
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LINEAR  AND  HOMOMORPHIC  FILTERING 

High  frequency  emphasis  filtering  of  the  FLIR  imagery  has  a two-fold  purpose 

• Attenuating  the  low  spatial  frequencies  reduces  the  global 
dynamic  range,  without  sacrificing  local  contrast 

• High  frequency  emphasis  crispens  edges  and  other  fine  detail 
that  might  not  have  been  otherwise  visible. 

The  above  concept  is  also  evident  in  the  local  area  gain  brightness  model, 
where  the  local  brightness  envelope  (see  Figure  15a)  denoted  by  the  broken 
line,  comprises  the  low  frequency  components  upon  which  the  higher  fre- 
quency local  variation  is  superimposed.  Attenuating  the  low  frequencies 
compresses  the  low  frequency  envelope  thereby  reducing  the  overall  dynamic 
range.  Of  course,  removal  of  the  low  frequency  components  is  not  desired. 
The  corresponding  filter  structure  and  the  overall  frequency  response  is 
shown  in  Figure  16a  and  b.  The  low  frequency  gain  A^,  the  high  frequency 
gain  A^,  and  the  3dB  frequency  fg^^g  of  the  low  pass  filter  should  be  tuned 
to  the  imagery. 

Homomorphic  filtering  is  implied  in  the  LOG  and  EXP  boxes.  Homomorphic 
processing  assumes  a multiplicative  model  of  the  image  (reference  3),  which 
is  true  of  imagery  in  the  visible  and  near  IR  parts  of  the  spectrum  where 
the  image  can  be  modeled  as  a product  of  low  frequency  illumination  and 
higher  frequency  reflective  components.  The  multiplicative  model  does 
not  seem  to  be  directly  applicable  for  FLIR  imagery,  but  can  be  easily 
incorporated  i/i  the  present  filter  structure. 
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The  Low  Pass  Filter 


The  basic  low  pass  filter  used  in  Figure  16a  can  be  realized  in  different 
ways.  We  are  designing  two  basic  realizations:  (1)  direct  convolution 
finite  impulse  response  (FIR)  Gaussian  filters  and  (2)  recursive  infinite 
impulse  response  (HR)  filters. 


Gaussian  Filter  (FIR) 


The  low  pass  filter  has  the  following  transfer  function: 


H(f^.fy)  = Exp 


■ .2  2 
f + f 
- _x y 

2a  ^ 


and  the  corresponding  FIR  convolution  filter  also  has  a Gaussian  shape: 

1 


h(x,y)  = 


2lf  r 


/2 


Exp 


2^2 
X + y 


2a 


where 


/ « 
a " 


271(7 


This  filter  does  not  have  a sharp  frequency  transition  but  has  the  advantage 
that  the  spatial  impulse  response  can  be  truncated  to  a finite  area,  without 
Gibb's  phenomenon  occuring.  Also,  the  filter  is  isotropic;  i.e.,  it  possesses 
circular  symmetry.  From  the  realization  standpoint,  it  has  the  advantage  of 
being  separable;  i.e,, 

h(x,y)  = h^(x)  hy(y) 
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Therefore,  a two-dimensional  convolution  can  be  cast  as  two  sets  of  succes- 
sive one-dimensional  convolutions; 


y<m,  n) 


M-l  M-1 
E E x(k.  1) 

k=0  1=0 

M-1  M-1 

E h (m-k)  E 

k=0  ^ 1=0 


h(m-k. 


x(k,  1) 


n-l) 


h (n-l) 

y 


where  m,  n = 0, 


N 


A direct  two-dimensional  convolution  of  a MxM  filter  with  a NxN  image 
2 2 

would  require  N M multiplications  and  additions.  This  separability 

enables  us  to  filter  with  one -dimensional  convolution  of  the  rows  of  the  image 

followed  by  one-dimensional  convolution  of  the  resultant  columns  as  implied 

2 

above.  This  requires  only  2N  M multiplications  and  additions--a  distinct 

2 2 

advantage  over  the  N M requirement,  especially  for  large  filters.  How- 
ever, the  filtered  rows  have  to  be  transposed,  so  that  the  columns  can  be 
filtered--a  disadvantage  when  working  with  a sequential  mass  medium  such 
as  magnetic  tape.  Therefore,  we  are  implementing  both  the  direct  2-D 
spatial  filter  for  small  filter  sizes,  and  the  separable  filter  for  wide  impulse 
response  filter  realization. 

Choice  of  3dB  Frequency 


The  choice  of  3dB  frequency  of  the  basic  low  pass  filter  in  Figure  16a  is 
very  important  in  contrast  enhancement.  If  only  edge  emphasis /crispening 
is  desired,  the  corresponding  filter  can  only  be  a few  pixels  wide.  The 
corresponding  frequency  response  is  shown  in  Figure  17a  where  the  3dB 
frequency  is  very  high.  On  the  other  hand,  for  attenuating  the  low  frequency 
components  for  dynamic  range  compression,  we  should  have  a much  lower 
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3dB  cutoff  in  the  filter,  Figure  17b,  which  requires  much  larger  spatial 
filters.  Unless  the  fast  Fourier  convolution  is  used,  direct  convolutional 
filtering  is  computationally  expensive  for  large  filters.  Hence,  the  separable 
Gaussian  filter  will  be  useful  for  this  implementation. 

Recursive  Filters 


The  large  impulse  responses  needed  for  spatial  filtering  for  contrast  enhance- 
ment can  be  realized  using  simple  two-dimensional  recursive  filters.  The 
advantages  of  recursive  filters  are: 

• Computational  efficiency--any  desired  3dB  frequency  (corre- 
sponding to  large  impulse  responses)  is  realizable  with  a low 
order  recursive  filter. 

• A two-dimensional  recursive  filter  can  be  designed  as  a product 
of  two  one -dimensional  low  pass  filters.  The  resultant  filter, 
although  not  isotropic,  may  be  quite  adequate  for  contrast 
enhancement  (reference  4), 

• Ease  of  implementation- -looking  ahead  to  2nd  Generation  FLIR 
implementation,  low  order  (1  or  2)  two-dimensional  recursive 
filters  are  readily  implementable  using  present  CCD  technology. 

Only  2n  line  buffers  are  needed  for  an  n*'  order  recursive 
filter.  It  is  expected  that  n»l  or  2 will  be  sufficient  for  the 
present  application. 

Low  order,  two-dimensional  recursive  Butterworth  filters  are  being  pro- 
grammed to  this  end.  If  these  filters  prove  to  be  of  value  for  contrast 
enhancement,  they  will  provide  ready  CCD  implementation. 
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Relation  Between  LAGBC  and  Linear  Filtering  for  Contrast  Enhancement 

The  basic  local  area  brightness  control  shifts  the  local  means  closer  to  the 
global  mean  to  reduce  the  overall  dynamic  range;  the  contrast  term  enhances 
the  local  contrast.  This  was  illustrated  in  Figure  15.  As  pointed  out  before, 
this  corresponds  to  attenuating  the  low  spatial  frequencies  and  enhancing  the 
high  spatial  frequencies.  This  can  be  seen  as  follows.  We  recognize  that 
the  local  mean  is  the  convolution  of  the  rectangular  window  with  the  image. 

In  the  spatial  frequency  domain,  this  corresponds  to  filtering  with  a low  pass 
filter  that  is  a product  of  two  Sine  functions;  i.e. , 

H(f  , f ) = Sinc(f  ) Sinc(f  ) 

X y X y 

We  may  as  well  replace  this  filter  by  a Gaussian  (or  other)  low  pass  filter 
with  the  same  noise  equivalent  bandwidth.  Figure  18  is  a realization  of  the 
LAGBC  using  such  a low  pass  filter  instead  of  the  averaging  process.  The 
frequency  response  at  each  stage  of  the  process  is  shown  in  the  insets.  For 
dynamic  range  compression,  the  gain  A is  less  than  unity  and  the  high  fre- 
quency gain  G is  greater  than  unity.  The  resultant  response  is  much  like 
the  basic  high  frequency  emphasis  filter  described  above.  The  spike  at  the 
origin  corresponds  to  the  bias  added  to  restore  the  average  brightness. 

There  are  differences;  for  example,  the  high  frequency  gain  G is  a function 
of  the  low  frequency  component  which  makes  the  process  non-linear.  Never- 
theless, the  remarkable  similarity  between  the  two  approaches  leads  us  to 
believe  that  they  are  equivalent  in  function.  This  exercise  serves  a two-fold 
function; 

• Relates  the  heuristic  LAGBC  concept  to  the  rigorous  linear 
filter  theory.  The  cutoff  frequency  of  the  low  pass  filter--the 
key  parameter- -determines  the  size  of  the  averaging  window 
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in  LAGBC  and  vice  versa.  The  low  frequency  gain  A is 
determined  by  the  range  of  the  scene,  and  the  high  frequency 
gain  G,  by  the  average  local  intensity. 

• Ease  of  implementation--computation  of  local  mean  in  LAGBC 
(especially  over  large  windows)  may  be  computationally  unfea- 
sible. In  the  equivalent  LPF  realization,  Figure  18,  we  can 
implement  a simple  two-dimensional  first  order  recursive 
filter  which  can  be  readily  implemented  with  present  CCD 
technology. 
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Figure  18.  A Linear  Filter  Approach 
to  Gain -Bright ness  Control 
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SECTION  III 
MRT  ENHANCEMENT 

There  are  two  approaches  to  improving  system  sensitivity  by  image  pro- 
cessing; intra-frame  and  inter-frame  spatial  averaging.  In  intra-frame, 
smoothing  one  trades  resolution  for  enhanced  S/N,  In  inter-frame 
smoothing  the  tradeoff  is  basically  between  data  rate  and  S/N,  although 
resolution  will  be  lost  if  registration  is  not  maintained  between  stacked 
image  frames.  To  date,  only  intra-frame  smoothing  has  been  addressed 
in  the  program.  Median,  hysteresis  and  variable  width  spatial  filtering  are 
currently  being  investigated. 

The  n*^-dian  filter  simply  computes  the  median  intensity  in  a small  sliding  . 

window  (3x3  or  5x5)  and  replaces  the  pixel  intensity  at  the  center  with  the  * 

median  value  for  the  window.  This  filter  has  the  advantage  that  median  t 

values  are  measured  values,  hence  original  data  is  almost  preserved, 
albeit  rearranged.  It  appears  to  be  most  useful  as  a pre-filter  operation 

to  remove  large  noise  spikes.  ] 

I 

i 

Hysteresis  filtering  is  a method  which  preserves  significant  extrema  in  an 
image  while  smoothing  out  local  fluctuations.  Two-dimensional  hysteresis 
filtering  is  achieved  simply  by  applying  one-dimensional  filters  sequentially 
in  orthogonal  directions  (separability  assumption).  For  one  dimension,  the 
algorithm  is: 


r 


If  (k=l),  = Xj^ 

*k+l  ■ \+l  *k+l  ' 

Vi ' Vi  * *''2 

Otherwise,  = y^ 

Here,  x.  is  the  original  pixel  intensity,  y,  is  the  filtered  intensity,  and  L 

K K 

is  the  hysteresis  window  width. 

A number  of  variable  width  spatial  filter  algorithms  have  been  proposed, 
differing  mainly  in  the  criteria  used  to  select  the  filter  size.  Examples  are 
criteria  based  upon  the  local  gradient;  or  the  texture,  to  separate  object 
boundaries  and  avoid  smoothing  across  them;  or  upon  the  local  average 
intensity  to  make  the  smoothing  window  S/N  dependent  for  quantum  noise 
limited  pictures. 

In  FLIR  imagery  processing,  object  edge  detection  is  the  most  attractive 
criteria  because  object  shape  and  contrast  appear  much  more  useful  than 
texture  as  classification  clues.  Also,  although  FLIR  detectors  are  generally 
BLIP,  the  actual  range  of  S/N  in  the  scene  does  not  vary  markedly.  Hence, 
local  S/N  variation  is  not  a good  criteria  for  chosing  local  filter  width.  The 
range  of  filter  width  selection  should  depend  upon  the  average  S/N  in  the 
scene,  however,  if  effective  smoothing  is  to  be  accomplished. 

Use  of  the  gradient  as  a criterion  for  selection  of  a variable  filter  width  to 
avoid  smoothing  across  edges  leaves  something  to  be  desired.  The  reason 
is  that  no  error  in  smoothing  an  edge  results  if  the  edge  gradient  is  constant 
within  the  width  of  the  smoothing  window. 
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To  investigate  this  further,  assume  that  an  edge  has  the  form  of  a plane 
ridge,  at  least  over  the  smoothing  window  area.  If  the  x-coordinate  is  in 
the  direction  of  the  gradient,  the  edge  can  be  described  by 

2 3 

f(x,  y)  = a + bx  + cx  + dx 
The  smoothed  intensity  (for  a box  filter)  is  then 

x+L/2  y+L/2  2 2 ^ 3 

F(x,  y)  =— *•  fdu  rdv  f(u,  v)  = (a  + ) + (b  + -j— )x  + cx  + dx 
L x-L/2  y-L/2 

For  specified  error  in  the  smoothed  value  at  the  center  of  the  window 


F(0,0)  - ((0.0) 

C 

f(0,0) 

5 « It  T2 

a 

Thus,  the  smoothing  error  does  not  depend  upon  the  local  gradient,  b,  but 
rather  on  the  local  curvature,  c.  Furthermore,  we  now  have  a criterion  for 
selecting  the  smoothing  window  size  to  achieve  specified  error. 

This  analysis  suggests  the  following  approach.  At  each  pixel,  the  least 
square  error  estimate  of  the  coefficients  a and  c in  f(x,y)  is  computed  under 
the  assumption  that  x is  in  one  of  the  four  principle  directions  (0°,  90^, 

± 45°).  For  a 7 pixel,  1-D  window  these  are 

a = (7gg  - gj.)/2l  c = (g^  - 4g^)/84 

where 


3 

= Z 
i = -3 


f. 


^c  = 


3 

T i 
i = l 


f J 


r = f(x+iAx,y) 
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r ^ 


I 

I 

! Use  of  a 7 pixel  window  provides  the  advantage  of  some  noise  smoothing  in 

the  filter  width  selection  criterion. 

t 

! 

Having  computed  a and  c for  each  of  the  4 principle  directions,  the  maximum 
L for  which  | c|l^  s I2e|a|  is  found  in  each  direction.  Then  the  minimum  of 
the  max  L for  the  filter  width  is  chosen.  A variation  on  this  scheme  is  to 
orient  an  asymmetric  filter  along  the  edge  by  allowing  the  filter  width  in  the 
direction  normal  to  that  of  the  min  max  L to  remain  the  max  L for  that  direc- 
tion. 

These  variable  width  filters  are  being  investigated  using  2-D  Gaussian  filters 
of  equivalent  rectangular  width  L=0,  2,  3,  5,  and  7.  These  correspond  to 
S/N  gains  against  white  noise  of  0,  3,  5,  7 and  8.5  dB,  respectively. 

Figure  19  shows  results  of  the  addition  of  median  and  hysteresis  filtering  of 
a FLIR  image  with  white  Gaussian  noise  at  an  average  S/N  = 5dB  (peak  signal 
to  rms  noise).  The  dark  lines  are  the  result  of  a software  anomaly  which 
caused  the  Optronics  film  writer  occasionally  to  skip  lines.  In  Figure  19, 

A is  the  original  image,  B is  the  noisy  image,  C is  the  median  filtered 
image  (5x5  window)  and  D is  the  hysteresis  filtered  image  (L  = 8).  Both 
filters  preserve  target  edges,  but  do  not  reduce  the  grainy  quality  of  the 
picture  appreciably.  Computation  of  the  mean  square  error  between  the 
original  and  the  noisy,  median  and  hysteresis  filtered  images,  respectively, 
yields  an  overall  S/N  gain  of  about  4 dB  for  the  median  filter,  and  less  than 
1 dB  for  the  hysteresis  filter.  These  measures  correlate  with  the  differences 
in  subjective  quality  of  the  two  filtered  images. 


I 
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SECTION  IV 


RESOLUTION  RESTORATION 


Work  on  full  frame  focus  restoration  and  super  resolution  is  in  its  prelim- 
inary stages.  Space  variant  blur  functions  have  been  modeled  and  images 
blurred  for  subsequent  restorations.  However,  no  2-D  restoration  algorithms 
are  operational  at  this  time. 

The  question  of  depth  of  field  has  a bearing  on  the  OTF  model  for  resolution 
restoration  studies.  In  particular,  if  the  autofocus  subsystem  can  achieve 
focus  at  the  hyperfocal  distance,  depth  of  field  problems  are  greatly  alle- 
viated. 

Exact  geometrical  optics  analysis  gives  the  maximum  range,  r^  and  mini- 
mum range,  for  which  the  diameter  of  the  angular  blur  circle  is  w as: 

" l-(wr7D)  ^2  = 1-Kwr/D) 

where  r is  the  range  at  which  the  system  is  in  focus  and  D is  the  aperture 
diameter.  At  the  hyperfocal  distance 

r = D/w 
o 

the  depth  of  field  extends  over  r^/2  i r s ® . For  example,  with  a NFOV 
8-inch  aperture  and  w=0.05  mrad.  the  depth  of  field  extends  from  2km  out 
to  infinity.  With  a WFOV  2-inch  aperture  and  w=0.2  mrad,  r^/2  = l27m. 

This  implies  that  only  the  relatively  slight  space  variant  off-axis  blur  of  the 
OTF  need  be  restored  if  system  focus  at  the  hyperfocal  distance  is  main- 
tained. 
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In  case  r^r  , the  depth  of  field  6 is 
o 

, 2wr^/D  _ 

6 - •'j  ~ ^2  2 " ^ 

l-(wr/D)  r«D/w 


2 2 2 # 
2wr  2w  r 


for  a detector  of  width  "a"  and  IFOV=w.  Thus  it  is  only  at  very  short  ranges 
where  the  depth  of  field  is  severely  constrained  by  F#  and  resolution. 


SECTION  V 


STATISTICAL  ANALYSIS  OF  IMAGERY 


Statistical  analysis  of  FLIR  imagery  is  being  performed  to  better  quantify 
the  characteristics  of  thermal  imagery,  and  to  provide  quantifiable  measures 
of  the  quality  of  enhancement  processes.  The  statistics  being  extracted  can 
be  categorized  as  shape,  intensity  and  texture  features.  Each  is  discussed 
separately  below. 

SHAPE  STATISTICS 


All  shape  statistics  are  derived  from  an  object-boundary  tracing  routine  and 
the  resulting  chain  code  which  describes  the  perimeter  of  the  object.  The 
boundary  extraction  routine  uses  the  Sobel  gradient  operator.  Referring  to 
Figure  20,  the  Sobel  gradient  operator  is; 

G (E)  = S + S 
s X y 

where 

s = |aa-ab|  s = |al-au| 

X ' y 

AA  = C + 2F  + I,  AB  = A + 2D  + G ■ 

AL  = G + 2H  + I,  AU  = A + 2B  -I-  C 


A 

B 

C 

D 

E 

F 

G 

H 

I 

Figure  20.  Sobel  Gradient  Geometry 
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To  find  a first  point  on  the  boundary,  an  initial  point  (Figure  21)  on  the 
outside  of  the  target  is  specified,  along  with  the  direction  D^,  to  proceed  so 
as  to  intercept  the  boundary.  At  a point  the  Sobel  gradient  exceeds  a pre- 
determined threshold.  Rather  than  start  tracking  the  boundary  at  this  point, 
however,  two  new  paths,  and  parallel  to  and  stradling  Pq,  are  chosen 

and  points  Q and  Q„  are  found  on  the  boundary  similar  to  Qq. 

1 2 


Figure  21.  Initial  Boundary  Search 

From  Qq,  and  Q2  and  their  eight  neighbors,  we  compute  an  average 
boundary  intensity  and  an  average  boundary  Sobel  gradient.  Qq  is  chosen 
as  the  starting  point  for  tracing  the  boundary  in  a clockwise  direction.  The 
next  boundary  point  is  determined  from  the  values  of  G,  Sx,  Sy,  AA,  AB, 

AL  and  ALT  of  that  point,  as  well  as  the  average  boundary  intensity  and  Sobel 
gradient  following  the  logic  illustrated  in  Figure  22, 

This  procedure  does  not  necessarily  terminate  at  the  starting  point  Qq. 
Typically  the  boundary  closes  as  below.  By  discarding  the  leading  edge 

^(^0' ® 


closed  boundary  starting  and  ending  at  Rq  is  generated. 
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Figure  22.  Boundary  Tracing  Logic 


1 


Figure  23  shows  a line  printer  plot  of  a tank  extracted  by  this  technique. 
Pixel  intensities  interior  to  the  target  are  negated  for  ease  of  separating 
target  from  background  when  extracting  brightness  and  texture  features. 

The  notation  (10x20/137)  gives  the  (Height  x Width/Area)  of  the  target  in 
pixels . 

Figure  24  shows  an  example  of  the  difference  in  the  extracted  boundary  be- 
tween unfiltered  and  a median  filtered  image.  Median  filtering  helps  to 
eliminate  edg_  ctnomalies  which  can  on  occasion  cause  the  boundary  detector 
to  trace  the  edge  one  pixel  into  the  interior. 
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Figure  25  delineates  the  types  of  statistics  which  can  be  extracted  from  the 
boundary  chain  code.  The  chain  code  is  edited  to  smooth  rough  boundaries, 
and  partitioned  into  edges  based  upon  an  average  straightness  criterion. 

Size  independent  statistics  such  as  the  ratio  of  perimeter  to  root-area  (P/^), 
number  of  edges,  distribution  of  normalized  edge  lengths  (length/perimeter), 
and  differential  slope  (change  in  slope  from  edge  to  edge  as  the  boundary  is 
traced  in  a clockwise  direction)  are  most  significant  for  comparisons  amongst 
target  classes  over  an  ensemble  of  images.  Note  that  P/^A  ^ 4 for  the  6 
tanks  is  typical  of  square  objects,  as  opposed  to  more  or  less  circular  blobs 
where  P/^  ~|  2/^  would  be  expected. 

Figures  26  and  27  are  sample  histograms  of  average  intensity  difference 
(contrast)  and  gradient  of  the  boundary  obtained  from  the  boundary  extraction 
routine.  It  is  anticipated  that  these  two  features  may  provide  useful  infor- 
mation in  classifying  target  characteristics  in  the  future. 


BRIGHTNESS  STATISTICS 


With  the  target  separated  from  the  background,  intensity  histograms  for  the 
target  (Figure  28).  and  local  background  (Figure  29)  are  easily  extracted. 
From  these,  histogram  moments  and  contrast  measures  such  as: 


Average  Contrast: 
and 

Peak  Contrast: 


|M.p-M  1 

Q = 

A max(M,p,  M^) 

^ ^ ^T,  max'^B 

P max(I, 


T,  max'  ^B^ 


can  be  computed.  Here,  T refers  to  the  target,  B to  the  background,  M is 
a mean,  and  I—  is  the  peak  target  intensity. 
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DIFFERENTIAL  SLOPE 

28.6 
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Figure  30. 


Mean  vs.  Standard  Deviation  of  Brightness  Histograms 
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TEXTURE  FEATURES 


These  are  measures  to  quantify  the  grey  level  variations  in  the  image,  called 
texture.  We  have  implemented  two  classes  of  texture  definitions  to  evaluate 
FLIR  imagery; 

• Grey  level  difference  statistics  [Rosenfeld  (reference  1)]  based 
on  the  Haralick  (reference  2)  measured. 

• Texture  features  based  on  Max-Min  measures  developed  by 
Professor  O.  R.  Mitchell  of  Purdue  University. 

Grey  Level  Difference  Statistics 


Assume  that  the  texture  is  to  be  measured  over  a local  area  of  the  image. 
Consider  all  pairs  of  points  in  the  region,  exactly  at  a vector  distance 
6 = (Ax,  Ay)  apart  (as  in  Figure  31).  Let  |p(x+Ax,  y+Ay)  - p(x,y)|be  the 
grey  level  difference  for  each  such  pair.  The  histogram  of  these  grey  level 
differences  (of  all  pairs  of  points  exactly  6^  apart)  is  called  the  grey  level 
difference  histogram  p.(  ).  If  there  are  N grey  levels,  then  there  are  N 

g 

bins  in  this  difference  histogram.  This  is  a measure  of  the  probability 
density  of  grey  level  differences  occurring  in  the  image  at  a given  spacing 
and  orientation.  The  shape  of  this  histogram  is  a measure  of  the  texture. 
The  various  descriptors  (describing  this  histogram)  are: 


a. 

(Mean) 

N-1  2 

b.  r i p.(i) 
i»0  ^ 

(Contrast) 
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(x+Ax,  y+Ay) 


/ 6 = (ax, Ay) 


(x,y) 


5 = {(0.1),  (1.0),  (1,1),  (1,-1)} 
{(0.2).  (2.0),  (2.2),  (2,-2)} 


Figure  31.  Illustrating  Point  Pairs  for  Grey  Level  Difference  Statistics 


N-1 

c.  E Pj  (i)  (Angular  Second  Moment) 

i=0 

N-1 

d.  -E  p,(i)  Ln[p.(i)]  (Entropy) 

i=0  ® ® 

Figure  32  gives  the  grey  level  difference  histograms  for  a target  and  its 
immediate  background  area  for  6 = (0, 1),  (2,  2),  (4,  0)  and  (8,  8). 

The  difference  histograms  offer  a great  deal  of  information  on  the  texture. 
But  we  need  fewer  descriptors  of  texture  than  the  full  histograms.  Hence 
the  above  measures  (mean,  contrast,  etc. ) were  computed  from  these  histo- 
grams. These  were  computed  on  the  target  (tanks)  and  background  areas  for 
six  NVL  supplied  images;  for  6=  [(0,  A).  (A,0),  (A,  A),  (A, -A)]  corresponding 
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TARGET  BACKGROUND 


to  orientations  of  0°,  90°,  45°  and  -45°  and  spacings  of  A = 1,  2,  4 and  8. 

The  parameters  (mean,  contrast)  were  averaged  over  all  directions  for  a 

\ 

given  spacing,  A.  Figure  33  is  the  corresponding  plot  of  the  mean  measure 
for  the  various  values  of  A = 1,  2,  4 and  8 for  the  six  tank  targets  and  their 
respective  backgrounds.  Note  that  the  separation  in  this  feature  for  each 
target  and  its  background  increases  with  A as  might  be  expected.  But  un- 
fortunately the  measures  do  not  cluster  into  background  and  target  values. 
This  might  also  be  expected  because  the  background  can  possess  a variety 
of  textures.  However,  this  is  too  small  a sample  to  be  statistically  mean- 
ingful. Figure  34  is  a similar  plot  for  the  contrast  measure. 
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Figure  33.  The  Mean  Measure  for  Various  Separations  6 


65 


f 


O TARGET 
O BACKGROUND 


□ 


6 


Figure  34.  Contrast  Measures  for  Various  6 
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Max-Min  Texture  Descriptors 


r 


The  general  idea  here  is  that  the  number  of  extrema  (i.  e. , maxima  and 
minima)  per  unit  length  along  a given  direction  in  the  image  is  a measure  of 
the  image  texture.  Since  we  are  dealing  with  irregular  areas  in  the  image 
(such  as  target/background)  this  concept  can  be  generalized.  The  distances 
between  successive  extrema  are  measured  in  a given  direction  in  the  area 
under  consideration.  Figure  35  illustrates  this.  The  distances  are  measured 
and  histogrammed  along  the  rows  and  columns  separately.  The  average  dis- 
tance, the  variance  of  these  distances,  as  well  as  skewness  and  excess  of 
these  histograms  are  computed  from  the  distance  histograms.  The  rows 
and  columns  are  further  subjected  to  hysteresis  filtering  with  different  lags 
(see  the  section  on  MRT  enhancement);  the  extrema  are  computed  and  the 
distances  measured  for  each  hysteresis  lag.  Hysteresis  filtering  with  differ- 
ent lags  yields  extrema  of  different  magnitudes.  The  corresponding  dis- 
tances can  be  heuristically  related  to  spatial  frequency  concepts,  with  l/(mean 
distance)  corresponding  to  the  average  frequency. 

Figure  36  shows  the  histogrammed  distances  between  extrema  for  target  and 
background  areas  for  an  NVL  image  for  different  hysteresis  lags,  0,  4,  8 and 
12. 

The  mean  and  standard  deviation  of  these  distances  (measured  along  rows) 
for  each  of  11  targets  and  its  backgrounds  are  plotted  in  Figures  37  and  38  I 

for  lags  4 and  8 respectively.  Again,  note  that  the  target  and  backgrounds  ' 

do  not  cluster  in  the  texture  measures,  but  each  target  differs  in  texture 
from  its  background. 

f 
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Although  the  texture  measures  presented  above  do  not  appear  to  cluster  into 
background  and  target  values  directly,  they  do  yield  discriminatory  infor- 
mation between  the  target  and  its  background.  It  will  be  interesting  to  see 
how  these  texture  features  are  transformed  by  the  contrast  enhancement 
and  resolution  restoration  processes. 


Figure  35.  Illustrated  Distances  Between  Extrema 
Along  a Given  Direction  in  the  Image 
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Figure  36.  Histograms  of  Distances  Between  Extrema  for  a Target 
and  Its  Background,  for  Different  Hysteresis  Lags  0,  4, 
8 and  12 
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Figure  37. 


Figure  38. 
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Plot  of  Mean  and  Standard  Deviation  of  Distances 
Between  Extrema  for  Hysteresis  Lag  = 4 
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Plot  of  Mean  and  Standard  Deviation  of  Distances 
Between  Extrema  for  Hysteresis  Lag  = 8 
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